#include "Pch.h"
#include "Visualizations.h"
#include "Atom.h"
#include "Config.h"

ContextData context;

// Clears all data (except the context itself) without destroying it.
// For initialization and after destroying data.
void ContextData::Clear() {
	font = 0;
	for (int i=0; i<sizeof(MGL_SPHERE_LISTS)/sizeof(mglMesh*); i++) {
		MGL_SPHERE_LISTS[0][i] = 0;
	}
	MGL_CYL_MESH = 0;
}

void ContextData::Init(wxGLCanvas *canvas) {
	if (!font) {
		if (!wxglContext) {
			wxglContext = new wxGLContext(canvas);

			canvas->SetCurrent(*wxglContext);
	#ifdef GLEW
			// Should only be called once.  If called again, and fails on second call
			// after succeeding on first, would be in trouble.  Not likely even if
			// ends up being called again, anyways.
			glewInit();
	#endif
		}
		else {
			canvas->SetCurrent(*wxglContext);
		}
		mglMesh::Begin(config->drawMode, 0);
		MGL_CYL_MESH = mglMakeCylinderMesh(0.10, config->cylRes);
		mglMesh::End();
		font = new GLFont(wxEmptyString, 12, canvas);
	}
}

ContextData::ContextData() {
	wxglContext = 0;
	Clear();
}

// Destroys data if it exists.
void ContextData::CleanupSharedData(wxGLCanvas *canvas) {
	// Shouldn't be needed, but just in case.
	Init(canvas);
	canvas->SetCurrent(*wxglContext);
	int needClear = 0;
	if (font) {
		delete font;
		needClear = 1;
	}
	for (int i=0; i<sizeof(MGL_SPHERE_LISTS)/sizeof(GLuint); i++) {
		if (MGL_SPHERE_LISTS[0][i]) {
			MGL_SPHERE_LISTS[0][i]->Destroy();
			needClear = 1;
		}
	}
	if (MGL_CYL_MESH) {
		MGL_CYL_MESH->Destroy();
		needClear = 1;
	}
	if (needClear) Clear();
}

void ContextData::Destroy(wxGLCanvas *canvas) {
	if (wxglContext) {
		CleanupSharedData(canvas);
		delete wxglContext;
		wxglContext = 0;
	}
}

void ContextData::DrawSphere(int type, int id, Color c) {
	if (!MGL_SPHERE_LISTS[type][id]) {
		AtomicInfo *info = GetAtomicInfoById(id);

		if (!info) return;

		if (type == SUB_VAN_DER_WAALS_SPHERE) {
			MGL_SPHERE_LISTS[type][id] = mglMakeSphereMesh(info->vanDerWaalsRadius * 0.20, config->sphereRes);
		}
		else {
			MGL_SPHERE_LISTS[type][id] = mglMakeSphereMesh(info->vanDerWaalsRadius, config->sphereRes);
		}

		if (!MGL_SPHERE_LISTS[type][id]) return;
	}
	MGL_SPHERE_LISTS[type][id]->Draw(c);
}

void CleanupRenderChain(ChainRenderData *chain) {
	for (int r=0; r<chain->numResidues; r++) {
		if (chain->residues[r].cartoon) {
			chain->residues[r].cartoon->Destroy();
		}
	}
	chain->cartoonInitialized = 0;
	free(chain->residues);
	chain->residues = 0;
	free(chain->connections);
	chain->connections = 0;
}

void RenderData::Release() {
	refs--;
	if (refs) return;
	for (int c=0; c<numChains; c++) {
		CleanupRenderChain(chains+c);
	}
	free(this);
}

static void CalcCylFrom(Matrix *m, Vector *coords, Vector *coords2) {
	Vector d;
	subVect(&d, coords2, coords);

	Vector h = {0, 1, 0};
	//subVect(&d, coords2, coords);
	Matrix m1, m2, m3;

	// Adjust cylinder's height.
	m3 = identity;
	m3.M[1][1] = lengthVect(&d);

	// Adjust cylinder's angle.
	Quaternion q;
	normalizeVect(&d);
	double angle = acos(dotVect(&h, &d));
	mulVect(&q.v, normalizeVect(crossVect(&q.v, &d, &h)), sin(angle/2));
	q.w = cos(angle/2);
	quatToMat(&m2, &q);

	// Move cylinder.

	m1 = identity;
	m1.M[0][3] = coords->x;
	m1.M[1][3] = coords->y;
	m1.M[2][3] = coords->z;

	Matrix temp;
	mulMat(&temp, &m2, &m3);
	mulMat(m, &m1, &temp);
}

// Calculate connection cylinder info in advance, at the moment.
void InitializeChain(ChainRenderData *render, ManagedChain *chain) {
	// Already initialized.
	if (render->residues)
		return;
	render->numResidues = chain->numTotalResidues;
	render->residues = (ResidueRenderData*)calloc(chain->numTotalResidues, sizeof(ResidueRenderData));
	render->numConnections = chain->atomInfo[chain->numTotalAtoms].firstConnectionIndex;
	render->connections = (ConnectionRenderData*)malloc(sizeof(ConnectionRenderData) * render->numConnections);
	Matrix prevInv;
	for (int a=0; a<chain->numTotalAtoms; a++) {
		PDBAtom *atom = &chain->atoms[a];
		ExtendedAtomInfo *info = &chain->atomInfo[a];
		Vector *coords = &atom->pos;
		if (!a || atom[0].resIndex != atom[-1].resIndex)
			prevInv = identity;
		for (int i=info[0].firstConnectionIndex; i<info[1].firstConnectionIndex; i++) {
			Vector coords2 = chain->connectionInfo[i].destAtom->pos;
			addVect(&coords2, &coords2, coords);
			mulVect(&coords2, &coords2, 0.5);

			Matrix m;
			// Makes sure cylinders are rotated the same way.  As a result, triangles of each half of cylinders
			// meet exactly, modulo roundoff.
			if (chain->connectionInfo[i].destAtom < atom) {
				CalcCylFrom(&m, coords, &coords2);
			}
			else {
				CalcCylFrom(&m, &coords2, coords);
			}
			Matrix merged;
			mulMat(&merged, &prevInv, &m);
			invertMatFull(&prevInv, &m);
			mglConvertMatrix(&render->connections[i].cylMatrix, &merged);
		}
	}
}

RenderData *CreateRenderData(PDBManager *pdbManager) {
	RenderData *out = (RenderData*) calloc(1, sizeof(RenderData) + (pdbManager->numChains-1) * sizeof(ChainRenderData));
	out->numChains = pdbManager->numChains;
	out->pdbManager = pdbManager;
	out->refs = 1;
	out->Initialize();
	return out;
}

void RenderData::Uninitialize() {
	for (int c=0; c<numChains; c++) {
		CleanupRenderChain(chains+c);
	}
}

void RenderData::Initialize() {
	for (int c=0; c<numChains; c++) {
		InitializeChain(chains+c, pdbManager->chains[c]);
	}
}

void RenderData::DrawBallAndStick(int c, int r, unsigned int residueFlags, Color color) {
	ChainRenderData *chainRender = &chains[c];
	ManagedChain *chain = pdbManager->chains[c];
	int firstAtom = chain->residues[r].firstAtom;
	int lastAtom = chain->residues[r].lastAtom;
	if (firstAtom < 0)
		return;

	// All cylinders are the same color.  That is not the case for all atoms.
	if (residueFlags & RESIDUE_COLOR_ATOMS) {
		color = config->firstColor[ID_SECONDARY_COLOR_DEFAULT];
	}

	// Draw cylinders.  All same color, so more efficient to do them all first.
	int index = chain->atomInfo[firstAtom].firstConnectionIndex;
	glPushMatrix();
	for (int a = firstAtom; a <= lastAtom; a++) {
		// Don't draw backbone to backbone lines when cartoon mode also activated.
		if ((residueFlags & RESIDUE_RENDER_CARTOON) && (chain->atoms[a].flags & ATOM_BACKBONE)) {
			while (index < chain->atomInfo[a+1].firstConnectionIndex) {
				mglMultMatrix(&chainRender->connections[index].cylMatrix);
				if (chain->connectionInfo[index].destAtom->flags & ATOM_BACKBONE) {
					index++;
					continue;
				}
				context.MGL_CYL_MESH->Draw(color);
				index++;
			}
			continue;
		}

		while (index < chain->atomInfo[a+1].firstConnectionIndex) {
			mglMultMatrix(&chainRender->connections[index].cylMatrix);
			context.MGL_CYL_MESH->Draw(color);
			index++;
		}
	}
	glPopMatrix();

	// Draw spheres.
	glPushMatrix();
	Vector offset = {0,0,0};
	for (int a = firstAtom; a <= lastAtom; a++) {
		if ((residueFlags & RESIDUE_RENDER_CARTOON) && (chain->atoms[a].flags & ATOM_BACKBONE))
			continue;
		int id = chain->atomInfo[a].atomicIntSymbol;

		if (residueFlags & RESIDUE_COLOR_ATOMS) {
			color = GetAtomicColorByAtomIndex(id);
		}

		// Saving a push/pop each sphere like this doesn't seem to make much difference.
		// Does seem to make thing a couple percent faster with cylinders, but hardly
		// worth it, at least in my tests.
		Vector *v = &chain->atoms[a].pos;
		glTranslatef(v->x - offset.x, v->y - offset.y, v->z - offset.z);
		offset = *v;
		context.DrawSphere(SUB_VAN_DER_WAALS_SPHERE, id, color);
	}
	glPopMatrix();
}

void RenderData::DrawSpaceFill(int c, int r, unsigned int residueFlags, Color color) {
	ManagedChain *chain = pdbManager->chains[c];
	int firstAtom = chain->residues[r].firstAtom;
	int lastAtom = chain->residues[r].lastAtom;
	if (firstAtom < 0) lastAtom = -100;

	glPushMatrix();
	Vector offset = {0,0,0};
	for (int a=firstAtom; a <= lastAtom; a++) {
		int id = chain->atomInfo[a].atomicIntSymbol;

		if (residueFlags & RESIDUE_COLOR_ATOMS) {
			color = GetAtomicColorByAtomIndex(id);
		}

		Vector *v = &chain->atoms[a].pos;
		glTranslatef(v->x - offset.x, v->y - offset.y, v->z - offset.z);
		offset = *v;
		context.DrawSphere(VAN_DER_WAALS_SPHERE, id, color);
	}
	glPopMatrix();
}

class MultiStripMeshMaker {
	int sides, shareVertices, numRings;
	int vertsPerRing;
	mglNormalVertex *ringVertices;

	int currentRingVerts;
	Vector *ring;
	Vector *ringNormals;
public:
	// sides must be even.
	MultiStripMeshMaker(int sides, int shareVertices) {
		this->sides = sides;
		this->shareVertices = shareVertices;
		numRings = 0;
		ringVertices = 0;
		vertsPerRing = sides * (1 + !shareVertices);

		ring = (Vector*) malloc(sizeof(Vector) * sides * 2);
		ringNormals = (Vector*) malloc(sizeof(Vector) * sides * 2);
		currentRingVerts = 0;
	}

	void Clear() {
		numRings = 0;
		currentRingVerts = 0;
	}

	~MultiStripMeshMaker() {
		free(ringVertices);
		free(ring);
		free(ringNormals);
	}

	void AddVertex(const Vector &vertex, const Vector &normal) {
		ring[currentRingVerts] = vertex;
		ringNormals[currentRingVerts] = normal;
		currentRingVerts++;
		if (currentRingVerts == 2*sides) {
			currentRingVerts = 0;
			AddRing(ring, ringNormals);
		}
	}

	void AddRing(const Vector *vertices, const Vector *normals) {
		if (numRings % 10 == 0) {
			// Last part of size is an overestimate of what is needed for the two caps.
			ringVertices = (mglNormalVertex*) realloc(ringVertices, sizeof(mglNormalVertex)*((numRings+10) * vertsPerRing + sides*2 + 10));
		}
		if (shareVertices) {
			for (int i=0; i<vertsPerRing; i++) {
				Vector normal = normals[2*i];
				addVect(&normal, &normal, &normals[(2*i + 2*sides - 1) % (2*sides)]);
				normalizeVect(&normal);
				for (int j=0; j<3; j++) {
					ringVertices[numRings*vertsPerRing+i].v3f[j] = vertices[2*i].v[j];
					ringVertices[numRings*vertsPerRing+i].n3f[j] = normal.v[j];
				}
			}
		}
		else {
			for (int i=0; i<vertsPerRing; i++) {
				Vector n = normals[i];
				normalizeVect(&n);
				for (int j=0; j<3; j++) {
					ringVertices[numRings*vertsPerRing+i].v3f[j] = vertices[i].v[j];
					ringVertices[numRings*vertsPerRing+i].n3f[j] = n.v[j];
				}
			}
		}
		numRings++;
	}
	mglMesh *MakeMesh() {
		if (numRings <= 1) return 0;
		GLushort *indices = (GLushort*) malloc(sizeof(GLushort) * (2 * (sides+2) + numRings * (2 * sides + 4)));
		int numVerts = vertsPerRing * numRings;
		int index = 0;
		// Currently doesn't handle odd numbers of sides.
		for (int i=sides/2-1; i>=0; i--) {
			int srcIndex[2];
			if (shareVertices) {
				srcIndex[0] = (vertsPerRing-i) % vertsPerRing;
				srcIndex[1] = 1*i+1;
			}
			else {
				srcIndex[0] = (vertsPerRing-2*i) % vertsPerRing;
				srcIndex[1] = 2*i+2;
			}
			for (int j=0; j<2; j++) {
				mglNormalVertex *v1, *v2;
				v1 = &ringVertices[numVerts];
				v2 = &ringVertices[srcIndex[j] + vertsPerRing];
				*v1 = ringVertices[srcIndex[j]];
				double sum = 0;
				for (int k=0; k<3; k++) {
					v1->n3f[k] = v1->v3f[k] - v2->v3f[k];
					sum += v1->n3f[k]*v1->n3f[k];
				}
				sum = 1/sqrt(sum);
				for (int k=0; k<3; k++) {
					v1->n3f[k] *= sum;
				}
				indices[index++] = numVerts++;
			}
		}
		int direction = vertsPerRing;
		for (int i=0; i<sides; i++) {
			int side1Base;
			int side2Base;
			if (shareVertices) {
				side1Base = i;
			}
			else {
				side1Base = i * 2;
			}
			side2Base = (side1Base + 1) % (vertsPerRing);
			if (direction < 0) {
				side1Base += (numRings-1) * vertsPerRing;
				side2Base += (numRings-1) * vertsPerRing;

				side1Base ^= side2Base ^= side1Base ^= side2Base;
			}
			if (i) {
				int lastIndex = indices[index-1];
				// Add degenerate triangle.
				indices[index++] = lastIndex;
				indices[index++] = lastIndex;
				indices[index++] = side1Base;
				indices[index++] = side1Base;
			}

			for (int j=0; j<numRings; j++) {
				indices[index++] = side1Base;
				indices[index++] = side2Base;
				side1Base += direction;
				side2Base += direction;
			}
			direction = -direction;
		}

		int lastIndex = indices[index-1];
		// Add degenerate triangle.
		// Could get rid of if drew after first side,
		// but then would have to completely change draw order
		indices[index++] = lastIndex;
		indices[index++] = lastIndex;
		indices[index++] = numVerts;
		indices[index++] = numVerts;
		// Currently doesn't handle odd numbers of sides.
		for (int i=0; i<sides/2; i++) {
			int srcIndex[2];
			if (shareVertices) {
				srcIndex[0] = vertsPerRing * (numRings-1) + (vertsPerRing-i) % vertsPerRing;
				srcIndex[1] = vertsPerRing * (numRings-1) + 1*i+1;
			}
			else {
				srcIndex[0] = vertsPerRing * (numRings-1) + (vertsPerRing-2*i) % vertsPerRing;
				srcIndex[1] = vertsPerRing * (numRings-1) + 2*i+2;
			}
			for (int j=0; j<2; j++) {
				mglNormalVertex *v1, *v2;
				v1 = &ringVertices[numVerts];
				v2 = &ringVertices[srcIndex[j] - vertsPerRing];
				*v1 = ringVertices[srcIndex[j]];
				double sum = 0;
				for (int k=0; k<3; k++) {
					v1->n3f[k] = v1->v3f[k] - v2->v3f[k];
					sum += v1->n3f[k]*v1->n3f[k];
				}
				sum = 1/sqrt(sum);
				for (int k=0; k<3; k++) {
					v1->n3f[k] *= sum;
				}
				indices[index++] = numVerts++;
			}
		}
		mglMesh *out = CreateMesh(ringVertices, numVerts, indices, index, GL_QUAD_STRIP, config->drawMode);
		free(indices);
		return out;
	}
};

void CreateCartoonMesh(ManagedChain *chain, ResidueRenderData *render) {
	Vector *temp = (Vector*) calloc(1, sizeof(Vector) * chain->numChainResidues*4 + sizeof(char) * (chain->numChainResidues+2)*2);
	Vector *points = temp;
	Vector *tangents = temp+chain->numChainResidues;
	Vector *uptangents = tangents+chain->numChainResidues;
	Vector *ups = tangents+chain->numChainResidues*2;

	// 1 extra on either end, so don't have to check ends in some cases.
	char *happyPrev = ((char*) (ups+chain->numChainResidues))+1;
	char *happy = happyPrev+chain->numChainResidues+2;
	for (int r=0; r < chain->numChainResidues; r++) {
		chain->residueInfo[r].flags &= ~(RESIDUE_HAS_SPLINE_AFTER|RESIDUE_HAS_SPLINE_BEFORE);
		PDBAtom *a = 0;
		double maxDist2 = 0;
		if (chain->residueInfo[r].flags & RESIDUE_IS_RESIDUE) {
			a = GetAtom(chain, r, " CA ");
			maxDist2 = MAX_ALPHA_CARBON_DISTANCE*MAX_ALPHA_CARBON_DISTANCE;
		}
		else if (chain->residueInfo[r].flags & RESIDUE_IS_DNA) {
			a = GetAtom(chain, r, " C4'");
			PDBAtom *a2 = GetAtom(chain, r, " C4'");
			maxDist2 = MAX_DNA_DISTANCE*MAX_DNA_DISTANCE;
		}
		if (!a) continue;
		happy[r] = 1;
		points[r] = a->pos;
		if (happy[r-1]) {
			Vector v;
			subVect(&v, &points[r], &points[r-1]);
			double dist2 = lengthSquaredVect(&v);

			if (dist2 > maxDist2) {
				continue;
			}

			// Don't let DNA/RNA match with AA
			if ((chain->residueInfo[r].flags^chain->residueInfo[r-1].flags) & RESIDUE_TYPE) {
				continue;
			}

			happyPrev[r] = 1;
		}
	}

	for (int r=0; r < chain->numChainResidues; r++) {
		if (!happy[r]) continue;
		if (!happyPrev[r] && !happyPrev[r+1]) {
			happy[r] = 0;
			continue;
		}

		if (happyPrev[r] && happyPrev[r+1]) {
			Vector d1, d2;
			subVect(&d1, &points[r-1], &points[r]);
			subVect(&d2, &points[r], &points[r+1]);
			mulVect(&tangents[r], subVect(&tangents[r], &points[r+1], &points[r-1]), 0.5);
			crossVect(&ups[r], &d1, &d2);
			if (lengthSquaredVect(&ups[r]) > 0.0000001) {
				normalizeVect(&ups[r]);
			}
		}
		else if (happyPrev[r]) {
			subVect(&tangents[r], &points[r], &points[r-1]);
		}
		else {
			subVect(&tangents[r], &points[r+1], &points[r]);
		}
	}

	for (int r=0; r < chain->numChainResidues; r++) {
		if (!happy[r]) continue;
		if (lengthSquaredVect(&ups[r]) < 0.0000001) {
			if (happyPrev[r]) {
				ups[r] = ups[r-1];
			}
			if (happyPrev[r+1]) {
				addVect(&ups[r], &ups[r], &ups[r+1]);
			}
			if (lengthSquaredVect(&ups[r]) < 0.0000001) {
				// Should only happen on 2 point lines.
				if (fabs(tangents[r].x) < 0.8) {
					ups[r].x = 4;
				}
				else {
					ups[r].y = 4;
				}
			}
		}
		normalizeVect(&ups[r]);
		if (happyPrev[r] && dotVect(&ups[r], &ups[r-1]) < 0)
			mulVect(&ups[r], &ups[r], -1);
		if (happyPrev[r] && chain->residueInfo[r].flags & RESIDUE_IS_ALPHA  && !(chain->residueInfo[r+1].flags & RESIDUE_IS_ALPHA)) {
			//mulVect(&ups[r], &ups[r], -1);
		}


		if (happyPrev[r] && happyPrev[r+1]) {
			mulVect(&uptangents[r], subVect(&uptangents[r], &ups[r+1], &ups[r-1]), 0.5);
		}
	}

	Vector *mpoints = (Vector*) malloc(sizeof(Vector) * 3 * (config->cartoonRes*2-1));
	Vector *mhoriz = mpoints + config->cartoonRes*2-1;
	Vector *mvert = mhoriz + config->cartoonRes*2-1;

	MultiStripMeshMaker *currentMesh = new MultiStripMeshMaker(4, 0);
	MultiStripMeshMaker *nextMesh = new MultiStripMeshMaker(4, 0);

	for (int r=1; r<chain->numChainResidues+2; r++) {
		if (r >= 2) {
			/*if (chain->residueInfo[r-2].flags & RESIDUE_HAS_SPLINE_BEFORE) {
				if (chain->residueInfo[r-2].flags & RESIDUE_HAS_SPLINE_AFTER) {
					render[r-2].cartoon = CreateMesh(currentSpline, 8*config->cartoonRes, doubleSplineIndices, numDoubleSplineIndices, GL_QUAD_STRIP, config->drawMode);
				}
				else {
					render[r-2].cartoon = CreateMesh(currentSpline, 4*config->cartoonRes, splineIndices, numSplineIndices, GL_QUAD_STRIP, config->drawMode);
				}
			}
			else if (chain->residueInfo[r-2].flags & RESIDUE_HAS_SPLINE_AFTER) {
				render[r-2].cartoon = CreateMesh(currentSplinePart2, 4*config->cartoonRes, splineIndices, numSplineIndices, GL_QUAD_STRIP, config->drawMode);
			}//*/
			render[r-2].cartoon = currentMesh->MakeMesh();
			MultiStripMeshMaker *tempMesh = currentMesh;
			currentMesh = nextMesh;
			nextMesh = tempMesh;
			nextMesh->Clear();
		}
		if (r > chain->numChainResidues)
			continue;
		if (!happyPrev[r]) continue;
		chain->residueInfo[r].flags |= RESIDUE_HAS_SPLINE_BEFORE;
		chain->residueInfo[r-1].flags |= RESIDUE_HAS_SPLINE_AFTER;

		double h[2] = {0.20, 0.20};
		double width = 0.20;

		if (chain->residueInfo[r].flags & RESIDUE_IS_DNA) {
			width = 0.5;
			h[0] = h[1] = 0.5;
		}
		else {
			for (int q=0; q<2; q++) {
				if (chain->residueInfo[r-1+q].flags & RESIDUE_IS_BETA) {
					h[q] = 1.5;
				}
				else if (chain->residueInfo[r-1+q].flags & RESIDUE_IS_ALPHA) {
					h[q] = 1.5;
				}
			}
		}

		for (int p=0; p<config->cartoonRes*2-1; p++) {
			double s = p/(double)(config->cartoonRes*2-2);
			double s2 = s *s;
			double s3 = s2*s;

			double p1c =  2*s3 - 3*s2     + 1;
			double t1c =    s3 - 2*s2 + s;
			double p2c = -2*s3 + 3*s2;
			double t2c =    s3 -   s2;

			double p1c1 =  2*s2*3 - 3*s*2;
			double t1c1 =    s2*3 - 2*s*2 + 1;
			double p2c1 = -2*s2*3 + 3*s*2;
			double t2c1 =    s2*3 -   s*2;

			Vector v;
			Vector dv;
			Vector temp;
			Vector side, up;
			static Vector lastDv, lastSide;
			for (int i=0; i<3; i++) {
				v.v[i] = p1c * points[r-1].v[i] + p2c * points[r].v[i] + t1c * tangents[r-1].v[i] + t2c * tangents[r].v[i];
				dv.v[i] = p1c1 * points[r-1].v[i] + p2c1 * points[r].v[i] + t1c1 * tangents[r-1].v[i] + t2c1 * tangents[r].v[i];
				temp.v[i] = p1c * ups[r-1].v[i] + p2c * ups[r].v[i] + t1c * uptangents[r-1].v[i] + t2c * uptangents[r].v[i];
			}
			double height = p1c * h[0] + p2c * h[1];
			if (lengthSquaredVect(&dv) < 0.000000001) {
				dv = lastDv;
			}
			else {
				normalizeVect(&dv);
				lastDv = dv;
			}
			crossVect(&side, &dv, &temp);
			if (lengthSquaredVect(&side) < 0.000000001) {
				side = lastSide;
			}
			else {
				normalizeVect(&side);
				lastSide = side;
			}
			//crossVect(&up, &side, &dv);
			up = temp;
			normalizeVect(&up);

			mpoints[p] = v;
			mvert[p] = up;
			mhoriz[p] = side;
			mulVect(&up, &up, height);
			mulVect(&side, &side, width);
			{
				int w = p - config->cartoonRes+1;
				Vector cross[4];

				addVect(&cross[0], &v, &up);
				addVect(&cross[1], &cross[0], &side);
				subVect(&cross[0], &cross[0], &side);

				subVect(&cross[2], &v, &up);
				subVect(&cross[3], &cross[2], &side);
				addVect(&cross[2], &cross[2], &side);

				for (int i=0; i<4; i++) {
					Vector v[2], n[2];
					v[0] = cross[i];
					v[1] = cross[(i+1)%4];
					switch ((i+0)%4) {
						case 0:
							n[0] = up;
							break;
						case 1:
							n[0] = side;
							break;
						case 2:
							mulVect(n, &up, -1);
							break;
						case 3:
							mulVect(n, &side, -1);
							break;
					}
					// if (!(i&1)) mulVect(n, n, 0);
					n[1] = n[0];
					if (p < config->cartoonRes) {
						currentMesh->AddVertex(v[0], n[0]);
						currentMesh->AddVertex(v[1], n[1]);
					}
					if (w >= 0) {
						nextMesh->AddVertex(v[0], n[0]);
						nextMesh->AddVertex(v[1], n[1]);
					}
				}
			}
		}
	}
	delete currentMesh;
	delete nextMesh;
	free(mpoints);
	free(temp);
}

void RenderData::DrawCartoon(int c, int r, Color color) {
	ResidueRenderData *render = &chains[c].residues[r];
	if (!render->cartoon) {
		if (chains[c].cartoonInitialized)
			return;
		CreateCartoonMesh(pdbManager->chains[c], chains[c].residues);
		chains[c].cartoonInitialized = 1;
		if (!render->cartoon)
			return;
	}
	render->cartoon->Draw(color);
}
